import numpy as np
import matplotlib.pyplot as plt

#===================================================================

T_H1 = 1.57801e5


def hici_hui(T):
    """ Fit to hydrogen collisional ionization from
    http://adsabs.harvard.edu/abs/1997MNRAS.292...27H """

    lam = 2 * T_H1/T
    num = 21.11 * lam**(-1.089)
    den = 1 + (lam/0.354)**(0.874)
    rate = num / den**(1.101)    
    rate = rate * np.exp(-0.5*lam)    
    rate = rate * T**(-1.5)
    return rate


def hiirca_hui(T):
    """ Fit to hydrogen recombination (caseA) from
    http://adsabs.harvard.edu/abs/1997MNRAS.292...27H """
    
    lam = 2 * T_H1/T
    num = 1.269e-13 * lam**(1.503)
    den = 1 + (lam/0.522)**(0.470)
    rate = num / den**(1.923)
    return rate


def hiircb_hui(T):
    """ Fit to hydrogen recombination (caseB) from
    http://adsabs.harvard.edu/abs/1997MNRAS.292...27H """
    
    lam = 2 * T_H1/T
    num = 2.753e-14 * lam**(1.500)
    den = 1 + (lam/2.740)**(0.407)
    rate = num / den**(2.242)
    return rate



if __name__ == '__main__':

    print 'hello' 

    rc('xtick', labelsize=20)
    rc('ytick', labelsize=20)

    N = 100
    logT_min = 3.6
    logT_max = 5.4

    logT = np.linspace( logT_min, logT_max, N )
    T = 10**logT

    CI = hici_hui(T)
    RCa = hiirca_hui(T)
    RCb = hiircb_hui(T)

    fig = plt.figure( figsize=(10,10) )
    ax = fig.add_subplot(111)

    ax.plot( logT, np.log10(CI/RCa), color='blue', lw=2.0 )
    ax.plot( logT, np.log10(CI/RCb), color='red', lw=2.0 )

    ax.set_xlim( (logT_min, logT_max) )
    ax.set_ylim( (-6,6) )

    ax.set_xlabel( r'$\log ( T \, {\rm [K]} )$', fontsize=26 )
    ax.set_ylabel( r'$\log ( \gamma / \alpha )$', fontsize=26 )

    plt.tight_layout()

    plt.savefig( 'gamma_over_alpha.png', dpi=100 )


